Diagrammatic determinantal quantum Monte Carlo methods: 
Projective schemes and applications to the Hubbard-Holstein model 
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We extend the weak-coupling diagrammatic determinantal algorithm to projective schemes as well 
as to the inclusion of phonon degrees of freedom. The projective approach provides a very efflcient 
algorithm to access zero temperature properties. To implement phonons, we integrate them out in 
favor of a retarded density-density interaction and simulate the resulting purely electronic action 
with the weak-coupling diagrammatic determinantal algorithm. Both extensions are tested within 
the dynamical mean field approximation for the Hubbard and Hubbard-Holstein models. 
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I. INTRODUCTION 

Diagrammatic determinantal quantum Monte Carlo 
(DDQMC), be it the weak-coupling expansion or hy- 
bridization expansion approach, is emerging as the 
method of choice for impurity solvers In compari- 
son to the Hirsch-Fye approach ^ they are continuous 
time methods and thereby free of Trotter errors, more 
efficient, and more flexible. In this article we concentrate 
on the weak-coupling algorithm. After a short review 
of our implementation of the algorithm, we show how 
to generalize it to projective schemes as well as to the 
inclusion of phonon degrees of freedom. 

Projective schemes have already been implemented in 
the framework of the Hirsch-Fye algorithm and used in 
the context of dynamical mean field theories 0, S 01- 
Very similar ideas for the formulation of a projective 
DDQMC algorithm may be used and are reviewed in Sec. 
mil With the projective DDQMC, we can reproduce re- 
sults of [Hi at a fraction of the computational cost and 
access much lower projection parameters. 

Phonon degrees of freedom have very recently been 
implemented in the hybridization formulation of the 
DDQMC Q. Since the hybridization approach is based 
on the expansion in the hybridization, the inclusion 
of phonons relies on a Lang-Firsov transformation. 
In the weak coupling approach it is more convenient 
to integrate out the phonons in favor of a retarded 
interaction. The purely electronic model may then be 
solved efficiently within the weak coupling DDQMC. 
In Sec. IIVI we present some details of the algorithm 
and provide test simulations for the Hubbard-Holstein 



model in the dynamical mean-field theory (DMFT) 
approximation. 

II. THE DIAGRAMMATIC DETERMINANTAL 
METHOD FOR HUBBARD INTERACTIONS 



Here we will briefly review the diagrammatic determi- 
nantal method for the Hubbard model 

Hu = -J2 ("^^T - 1/2) {n^^l - 1/2) , 



Ho 

(1) 

where ni_cr = c\^Ci^cr and {ci,a) creates (annihilates) 
a fermion in a Wannier state centered around site i (j) 
and with z-component of spin a. 

As will become apparent in subsequent sections, we 
rewrite the Hubbard interaction as 



U 



J2 E «T - 1/2 - s5) «^ - 1/2 + s5) 



(2) 



i s=±l 



to avoid the negative sign problem at least for impurity 
and one-dimensional models. After carrying out the sum 
over the Ising spins, s, one recovers the original Hubbard 
interaction up to a constant. As will be seen below an 
adequate choice of S to avoid the sign problem for a one- 
dimensional chain reads S = ^ -I- 0"*". 

A weak coupling perturbation expansion yields for the 
partition function: 



dr„ ^ n(["-ii,<T(n) - aa{si)] ■ ■ ■ [rii„,o-(r„) - a<j(s„)])o ■ (3) 



Here, we have defined 

a^(s) = 1/2 asb 



and (•)n = Tr 



/Zq with Zq — Tr 



(4) 
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Note that the Ising field s has obtained an additional time 
index. The thermal expectation value is the sum over all 
diagrams, connected and disconnected, of a given order 

I 



with Green's functions 

G°,(ri,r2) = (rcl(ri)c,(r2))o, (6) 

which we have assumed to be spin independent. In the 
above, T corresponds to the time ordering. Defining a 
configuration, C„, by the n Hubbard vertices, as well as 
the Ising spins introduced in Eq. ([2]) 

C„ = {[ii, Ti, si] • • • [i„, r„, s„]} , (7) 

and the sum over the configuration space by 

E = E/ dnE--7 dr„^, (8) 

C„ 71=0-^^ ii.si "^0 i„,s„ 
I 



For any given configuration of vertices C„, Wick's theo- 
rem holds. Hence, any observable can be computed from 
the knowledge of the single particle Green's function 



71 

- E Gl^ir,n^){M-')^^^Gl^^iTs,n). (12) 

r,s—l 



Here, we have assumed that the non-interacting Green's 
functions are spin independent. As a consequence of 
the above equation, it becomes apparent that one can 
measure directly the Matsubara Green's functions. This 
aspect facilitates the implementation of the algorithm 
within the framework of dynamical mean-field theories. 



n. Using Wick's theorem this sum can be expressed as a 
determinant where the entries are the Green's functions 
of the non-interacting system. 



(5) 



I 

the partition function can conveniently be written as 

^=E(-|)"n^'^*^^'^(^")- 

Here is the n x n matrix of Eq. Obscrvables, 
0(t), can now be computed with 

^ Ec„ (-|)"n.detM.(C„)((6(r)))c„ 
^ ^"^^ Ec„(-f)"n.detA./.(aO 

where for 0(r) — Yl^ Oa{T) we have 



(11) 

I 

A. Sign problem 

In auxiliary field determinantal methods @ it is known 
that the presence of particle-hole symmetry can be used 
to avoid the negative sign problem. An identical state- 
ment holds for the diagrammatic determinantal method. 
We assume that Hq is invariant under the particle-hole 
transformation 

ia ^ (-I)'C,,.. (13) 

As a consequence, 

{TY\_n[ni^,a{Tr) " aa{sr)])o ^ (14) 

r=l 

(-l)"(r Yl n K,.(r,) - a-.(s,)])o 

r=l 



(T [h^j^{Ti) - acr(si)] • • • [n<T,,;„(T„) - ao-(s„)])o = 

/ G°,.»i(n,Tl) - a^isi) G0^,.,(Tl,T2) 



det 



V G?„,.,(r„,Ti) 



) \ 



UAT\ 



r(Ti) - aaisi)] ■ ■ ■ [ni„,o-(T„) - acr(s„)] da{T))o 



ncr(^["ii,'^('^i) ~ a<T(si)] • ■ ■ [ni„,o-(T„) - acr(Sn)])o 
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such that detM| = (-l)"detM|. A ghmpse at Eq. ^ 
wiU confirm the absence of sign problem for this special 
case. The above results is independent on the choice of 
S introduced in Eq. As we will see in Sec. Ill CI 

the algorithm is optimal at i5 = 0. In this special case, 
dctM| = detMx = (-l)"detMx such that only even 
values of n occur in the sampling. We note that this 
vanishing of the weight for odd values of n can be avoided 
by choosing a small value of 6. 

In one dimension and in the absence of frustrating 
interactions, there is no negative sign problem [l^. 
The diagrammatic approach also satisfies this property, 
provided that we choose (5 = 1/2 + 0"'". The quantity 
n^det M„(C„) in Eq. ^ is nothing but 

Tr[e-'^^« J][n,,,,(Ti)-a,(si)] 

••• [n,„,,(r„)-a,(s„)]]/Tr[e"'^^«] , (15) 

which we can compute within the real-space world-line 
approach [l^, Here, each world line configuration 

has a positive weight. Let us consider an arbitrary 
world-line configuration, and a site (i,T) in the space- 
time lattice. Irrespective if this site is empty, singly or 
doubly occupied the expectation value of the operator 
Off ["j.^lT) — acr{s)] will take a negative value. Recall 
that we have set 5 = 1/2-1-0"'". Hence, for each world 
line configuration the expectation value of the opera- 
tor [i'>T'i,criT) - acrisi)) ■ ■ ■ {ni„,a{Tn) ^ aa{Sn))] haS a 

sign equal to (—1)". Summation over all world line con- 
figurations yields the expression in Eq. (|15[) which in 
turn has a sign (—1)". This cancels the sign of the fac- 
tor (— [//2)" in Eq. ([5]), thus yielding an overall positive 
weight. 

In the rewriting of the Hubbard term (see Eq. ([2])) we 
have introduced a new dynamical Ising field so as to avoid 
the negative sign problem at least for the one-dimensional 
Hubbard model. Alternatively, one can choose a static 
Ising field and compensate for it by a redefinition of Hq. 
Such a static procedure is introduced in [l| . For the class 
of models considered, we have not noticed substantial 
differences in performance between static and dynamical 
choices of Ising fields. We however favor the dynamical 
version since it allows one to keep the SU(2) spin invari- 
ant form of the non- interacting Hamiltonian Hq. 



W{C) corresponds to the weight of the configuration. To 
add a vertex _^fj = which corresponds to the 
fact that one has to pick at random an imaginary time in 
the range [0, /3], a site i in the range 1 . . . (with A^ the 
number of sites) as well as an Ising spin. The proposal 
probability to remove a vertex ?c„^i^(7„ ~ TITT 'corre- 
sponds to the fact that one will choose at random one of 
the n + 1 vertices present in configuration C„+i, hence 



P< 



= mm — 



Uf3N n^detA/^(C„+i) 
(n + 1) n.detA/,(C„) ' 
jn + l) nadetAf^(C„) 

C//3A n.detM,(C„+i)^ 



(17) 



Apart from the above addition and removal of vertices, 
we have implemented moves which flip the Ising spins at 
constant order n as well as updates which move Hubbard 
vertices both in space and time. 



C. Tests 

The efficiency of the approach relics on the autocorre- 
lation time, which has to be analyzed on a case to case 
basis, as well as on the average expansion order param- 
eter. For a general interaction term Hi the average ex- 
pansion parameter is given by 

1 (-l)"n 
{n} = —y -. — / dri • • • / dT„ 



Z ^ n\ 







/ dn-- 


■ / dT,„ / 


h 


Jo Jo 



Z ^-^ ml 



X (Ti?i(Ti)---i?i(r™)i?i(T))o 



MHiir)) . 



(18) 



For the Hubbard model and replacing Hi by the form of 
Eq. ([2]) we obtain: 

(n) =-/3C/^[((n,,T-l/2)(r^u-l/2))-'52] . (19) 



B. Monte Carlo Sampling 

In principle two moves, the addition and removal of 
Hubbard vertices, arc sufficient [l^. In the Metropolis 
scheme, the acceptance ratio for a given move reads 



Pc^c 



rpO 



M{C') 



rro 



,W{C) 



,1 



(16) 



where T^,_,p corresponds to the probability of propos- 
ing a move from configuration C" to configuration C and 



Using the same techniques as in auxiliary field QMC 
methods the CPU costs for the calculation of the 
acceptance probability and the update for the addition 
or removal of a vertex scales as -n? . As apparent from 
Eq. (|19p a sweep consisting of updating all n vertices 
results in an effort of n? . Even though in this method 
M^^ is far better conditioned than in the classic deter- 
minantal methods [1, [^, [l2|, it has to be recalculated 
from scratch after several updates which involves an ef- 
fort of the order rv^ [l^]. Hence, because of these two 
limiting factors the CPU time scales as (/3 UN)^ which 
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timatc of the speedup reads (3200/400)^ = 512. Away 
from particle-hole symmetry, the speedup is less impres- 
sive since we have to set 6 = 1/2 + 0+ to avoid the neg- 
ative sign problem. We have confirmed the above state- 
ments for the Anderson impurity model, 



(20) 



+U ld\d^ - 1/2] ld[di - 1/2 



with Hubbard interaction U and hybridization V. Here 
cj, {ck,cr) creates (annihilates) a fermion in the conduc- 
tion band at momentum k with spin a and the disper- 
sion relation e{k) ~ —2t{cos{kx) + cos{ky)). The opera- 
tors dj^ and d„ create, annihilates an impurity electron, 
respectively. Our results including comparison with the 
Hirsch-Fye algorithm are presented in Fig. [TJ 

Finally let us note that applying the method to a one- 
dimensional Hubbard model of length TV, yields a very 
poor performance in comparison to standard finite tem- 
perature BSS auxiliary field algorithms [l3| since those 
methods scale as PUN^ [§]. 
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III. GENERALIZATION TO PROJECTIVE 
APPROACHES 



FIG. 1: (Color online) Green's function G(r) for the Ander- 
son model of Eq. (|20p . (a) Particle-hole symmetric point, 
fi — 0. To avoid the vanishing of the weight at odd values of 
n we have used 5 — 0.1 for this simulation, (b) Away from 
particle-hole symmetry jj, — —0.5. For this higher tempera- 
ture we provide a comparison with the Hirsch-Fye algorithm 
with Trotter step Art — 0.1. Note that in both cases we have 
used a 32 X 32 square lattice to generate the non-interacting 
Green's function G'^{t). (n) corresponds to the average order 
expansion parameter. 



is precisely the same scaling as in the Hirsch-Fye ap- 
proach. Apart from the absence of Trotter errors the 
advantage of the method lies in a large pre-factor. In the 
very special case of a single impurity, N = I, and for a 
particle- hole symmetric Hamiltonian Hq, the speedup is 
dramatic. Particle-hole symmetry allows one to set 6 — 
such that (n) < fiU/A. To obtain this upper bound we 
have set the double occupancy to zero. Hence a simu- 
lation at U/t = 4 and pt = 400 has a maximal average 
order parameter (n) = 400. In a Hirsch-Fye approach, 
one could opt for a Trotter step Art = 1/8 and hence 
3200 Trotter slices which determines the size of the ma- 
trices involved in the simulations. Hence an underes- 



Projective approaches rely on the filtering out of the 
ground state, \'^o) fi'om a trial wave function |^'t), which 
is required to be non-orthogonal to |5'o) 



(*T|e-Q-f^|*T) 



(21) 



For convenience and simplicity, we will assume that | 

is the ground state of Hq, such that for a given value of 6 

the right hand side of the above equation can be written 

as 



(^T|e-®^|*T> 



lim . (22) 

/3o^oo Xre-'^o^oe-^-f^ 



With the definition Z,, 



Tr 



^ — PqHq ^ — QH 



and 



Zpft — Tr 
yields 



-{l3o+e)Ho 



a weak coupling expansion 



dT„ ^ n(["ii,o-(Ti) - Q;cr(si)] • • • [rii„.cr(T„) - a„{sn)\)pfi , 



(23) 
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FIG. 2; (Color online) DMFT calculation of the double oc- 
cupancy for the half-filled Hubbard model as a function of i) 
temperature 1//3 and ii) projection parameter 1/0. In the 
projective approach, observables are computed in the imagi- 
nary time range [O/4,30/4]. This rather large measurement 
interval explains the slight discrepancy with the data of Ref. 

for the projective code and at finite values of Q. How- 
ever, extrapolation to Q —> cx) where the measurement range 
becomes irrelevant yields results consistent with Ref. 



where 



/p,o 



= Tr 



,-(/3o+e)Ho, 



/Zpfi. The similarity to 

the finite temperature algorithm is now apparent. We use 
Wick's theorem to express the expectation value on the 
right hand side of Eq. ([23|) in terms of the product of two 
determinants. Taking the limit — > oo we obtain 



(^'tIb^^^I^'t) = lim 



Zri 



/3o— »oo Zpfi 

u 



E 



[] dot M,,p(C„) .(24) 



Here, the sum runs over all configurations C„ as defined 
in Eq. and Eq. ©. Note that P in Eq. (0) has to 
be replaced by 9. The matrices M^^p have precisely the 
same form as the matrices Af^, but since we have taken 
the hmit /?□ — > oo, the thermal non- interacting Green's 
functions have to be replaced by the zero temperature 
ones: 



^0 



(Tl,r2) = {^t\Tc\{ti)cj{t2)\^t) 



(25) 



Hence, as in the Hirsch-Fye approach, the step from a 
finite temperature to zero-temperature code is very easy 
and essentially amounts in replacing the finite tempera- 
ture non-interacting Green's functions by the zero tem- 
perature ones. However, there is an important difference 
concerning measurements: measurements of observables 
which do not commute with the Hamiltonian have to be 
carried out in the middle of the imaginary time interval 
[0,0] to avoid boundary effects 0]. We have tested this 
approach by reproducing Fig. 1 of the article Q where 
the projective Hirsch-Fye algorithm was incorporated in 



the DMFT self-consistency cycle. Fig. [3] shows the re- 
sults for the half-filled Hubbard model at U/t = 4.8, den- 
sity of states N{uj) = -:^^^\JW'^ — and band- width 
W = \t. Excellent agreement with the former Hirsch-Fye 
based results Q both at finite temperatures and in the 
limit O ^ oo were obtained. However, the diagrammatic 
approach allows to access much larger projection param- 
eters and/or lower temperatures. We refer the reader to 
Ref. 01 for the implementation of the projective formal- 
ism in the self-consistency cycle. 



IV. APPLICATION TO THE 
HUBBARD-HOLSTEIN MODEL 



Weak coupling DDQMC allows a very simple inclusion 
of phonon degrees of freedom. The path we follow here 
is to integrate out the phonons in favor of a retarded 
interaction, and then solve the purely electronic model 
with the DDQMC approach. Starting from the Hubbard- 
Holstein model with Einstein phonons we show how to 
integrate out the phonons, describe some details of the 
algorithm and then present results within the DMFT ap- 
proximation. 



A. Integrating out the Phonons 

The Hubbard- Holstein Hamiltonian we consider reads 



(26) 



Here, hi = hi_„ and the last two terms correspond 
respectively to the electron-phonon coupling, and the 
phonon-energy. The Hamiltonian is written such that for 
a particle-hole symmetric band, half-filling corresponds 
to chemical potential /i = 0. Opting for fermion coherent 
states 



(27) 



Ci^a being a Grassmann variable, and a real space repre- 
sentation for the phonon coordinates 



(28) 



the path integral formulation of the partition function 
reads 



(29) 



with 



Su = 



+ ?75]KT(r)-l/2)(n,,i(r)-l/2) 



5e 



L 

Jo i 



A/9,2 (r) , fc 



2 2 
+ gg,;(T)(n,(T) - 1) . 

In Fourier space, 



(30) 



(31) 



where fim is a bosonic Matsubara frequency, the electron 
phonon part of the action reads 



^ / dr5]e-("--"'=^)(n,(r)-l).(32) 



Gaussian integration over the phonon degrees of freedom 
leads to a retarded density-density interaction: 



dr dr'Y., ,[n,{r)-l]D°{i-3,r-r')[n,(T')-l\ ^ ^33^ 



For Einstein phonons the phonon propagator is diagonal 
in real space. 



J, T-T')-5,,,|^P(T-r') with 



P(r) 



Wo 



2(1 - e-''3'^o) 



Hence the partition function of the Hubbard-Holstein 
model takes the form. 

Z = j [dcUc] X (35) 

^-(Su- dr Jl' dr^T,, ^n,{r)-l\D\^-j,T-r')[n,[r^)-l]) _ 

In the anti-adiabatic limit, Yvcaujo^oo P{t) = 6(t) such 
that the phonon interaction maps onto an attractive Hub- 
bard interaction of magnitude /k. We are now in a 
position to apply the DDQMC algorithm by expanding 
in both the retarded and Hubbard interactions. 



B. Formulation of DDQMC for the 
Hubbard-Holstein model 



To avoid the minus-sign problem at least for the one- 
dimensional chains we rewrite the phonon retarded inter- 



0.4 
0.35 
Q 0.3 
0.25 

0.^ 



(a) 

U/t=l, co,/t=0.2, Pt=40 




'0 



0.1 0.2 ^1^03 0.4 



g/(2McOo) 



12r 



(b) 

9- U/t=l,cOQ/t=0.2, pt=40 



^ 6 



3^ 



0.1 



><'0.S 



0.4- 



U/t=l, C0g/t=0.2, Pt=40 



0.5 




0.2 ,„0.3 0.4 0.5 




0.1 0.2 ,0.3 0.4 0.5 

g/(2McoQ) ' 



FIG. 3; (a) Double occupancy, (b) local charge susceptibil- 
ity and (c) local spin susceptibility for the Hubbard-Holstein 
model in the DMFT approximation. 



action as: 



Hp{t) 



X 



i,cr,a' s=±l 

[■ni,cr{T) - a+{s)] [ni^„>{T') ~ a+(s)] 



(36) 



For each phonon vertex, we have introduced an Ising vari- 
able: s. Summation over this Ising field reproduces. 
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up to a constant, the original interaction. Since the 
phonon term is attractive the adequate choice of signs 
is ct+{s) = 1/2 + sS, irrespective of the spin cr and cr'. A 
similar argument as presented in Sec. Ill Al then guaran- 
tees the absence of a sign problem for chains. Following 
Eq. ([2]) we rewrite Hubbard term as 

To proceed with a description of the implementation of 
the algorithm it is useful to define a general vertex 

V(t)={i,t,o,t' ,o' ,sM . (38) 

where h defines the type of vertex at hand, Hubbard 
(& = 0) or phonon (6 = 1). For this vertex we define a 
sum over the available phase phase 




a weight 

«^[nr)]=<^M|-4aP(r-r')|^, (40) 

as well as 

7?[y(T)] =,5m'5<,,|(5<,,,^5(t-t') 

K,T (^) - "+ (•?)] (^) - «- (•?)] + (41) 
<5b,i {ni.„{T) - q;+(s)] {n^^„,(T') - a+(s)] . 

With the above definitions, the partition function can 
now be written as 

V = E(-l)" / dri E MViin)] ■ ■■ dr„ 

J2 ^ \Vn{T„)] {TH [Vi{ti)] ■ ■ ■H [K(t„)])o . (42) 

V„(r„) 

As for the Hubbard model, a configuration consists of a 
set vertices C„ = {Vi(ti ),..., T^(r„)}. For a given con- 
figuration the thermal expectation value maps onto the 
product of two determinants in the spin up and spin down 
sectors. The Monte Carlo sampling follows precisely the 
scheme presented in Sec. IH Bl namely the addition and 
removal of vertices. 



C. Application to the Hubbard-Holstein model in 
the dynamical mean field approximation 

We have applied the above algorithm to the Hubbard- 
Holstein model within the dynamical mean-field the- 
ory approximation. Wc use a semicircular density 
of states, N{u!) = :^\/W^/4 - iv^ with band- width 
W = At. Throughout this section, we set U/t = 1, 
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FIG. 4: (Color online) The normalized histogram of the lo- 
cal density p of the Hubbard-Holstein model at half-filling 
with U /t — 1, ujo = 0.2t and /3t = 40, for several values of the 
electron-phonon coupling g. 



luq = 0.2t, fi ~ 0, and use the finite temperature imple- 
mentation of the algorithm at /3t = 40. The choice /i = 
corresponds to half-filling. Fig. [3]^a) plots the double 
occupancy, D = (ni_|ni.|) as a function of the electron- 
phonon coupling. To compare at best with the results 
of Ref. [ist we write the phonon coordinates in terms 
of bosonic operators, Q = ^2Muj ('^ ^ ^^)' ^"^^ V^o^ our 
results as a function of 

g^g/^2Mu;a . (43) 



Comparison of the results of Fig. [3I^a) with those 
of Ref. show excellent agreement and a critical 

electron-phonon coupling for the transition to the bipo- 
laronic insulator at cjc — 0.45f. 

We have equally computed the local spin and charge 
susceptibilities: 

Xc= dr([7i,,|(T) ± n,.|(T)] ± n,.J) . (44) 
s Jo 

As can be seen in Fig. ^h) in the vicinity the transition 
local charge fluctuations grow substantially and local spin 
fluctuations (see Fig. [Sjc)) are suppressed. The suppres- 
sion of xs signals singlet pairing of polarons. The chemi- 
cal potential /i = sets the average particle number p = 1 
but the particle number itself oscillates strongly between 
an empty or doubly occupied site. This can be seen by 
taking histograms as shown in Fig. 2) In the vicinity 
of the transition a two peak structure corresponding to 
a doubly occupied or empty site emerges. Close to the 
transition it becomes increasingly hard to guarantee a 
symmetric histogram - corresponding to the particle-hole 
symmetry of the model - and the simulation ultimately 
freezes in the doubly occupied or empty state. 
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V. CONCLUSIONS 

We have presented two extensions of the diagrammatic 
dcterminantal method: projective schemes as well as the 
inclusion of phonons. In both cases we have tested suc- 
cessfully the approach for the Hubbard as well as for the 
Hubbard-Holstein models in the dynamical mean-field 
approximation. The inclusion of phonons is not limited 
to Einstein modes and in principle any dispersion relation 
can be easily implemented. One of the strong points of 
the weak-coupling DDQMC can be extended very easily 
to larger clusters and used as a solver for cluster exten- 



sions of dynamical mean- field theories [14[. The crucial 
issue here is the severity of the sign problem as the clus- 
ter size increases. This is an issue which would need to 
be answered on a case to case basis. 
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